%% Assignment 8.1
% CIE4440

close all;clear all;clc;

%% 1.	First estimate h0, by plotting a regression line through the values of the lower discharges for the period 1995/96. What is h0?
load data

figure(1)
cax = gca;
set(cax,'NextPlot','add');
xlabel('Q [m3/s]');
ylabel('h [m]');

% data points scattered
scatter(data.Q,data.H)

% Regression line plotted
[r,a,h0] = regression(data.Q',data.H');
regx=0:max(data.Q);
regy=a*regx+h0;
plot(regx,regy);

fprintf('the h0, the offset is estimated on %f m\nThe correlation coefficient is %f\n',h0,r);
%% 2.	10-log(?) Hereafter plot Log (Q) against Log(h-h0). This results in a certain correlation coefficient. Try to optimize this value, that will give a new h0 value. What is your new h0?
xdata=log10(data.Q');
ydata=log10(data.H');

figure(2)
cax     = gca;
set(cax,'NextPlot','add');
xlabel('log10(Q m3/s)');
ylabel('log10(h m)');

% data points scattered
scatter(xdata,ydata);

% Regression line plotted
[r,a,h0_]   = regression(xdata,ydata);
regx        = linspace(0,max(xdata));
regy        = a*regx+h0_;
plot(regx,regy);

h0          = 10^h0_;

fprintf('The h0, the offset is estimated on %f m, which is a lot lower than the not-log\nThe correlation coefficient is %f, which is the same as before\n',h0,r);

%% 3.	Plot the discharge measurements + historic floods Log (Q) against Log (h-h0).
xdata=log10([data.Q;data.floodQ]');
ydata=log10([data.H;data.floodH]');

figure(3)
cax     = gca;
set(cax,'NextPlot','add');
xlabel('log10(Q m3/s)');
ylabel('log10(h m)');

% data points scattered
scatter(xdata,ydata);

% Regression line plotted
[r,a,h0_]   = regression(xdata,ydata);
regx        = linspace(0,max(xdata));
regy        = a*regx+h0_;
plot(regx,regy);

h0          = 10^h0_;

fprintf('the h0, the offset is estimated on %f m, which is a lot lower than the not-log\nThe correlation coefficient is %f, which is the same as before\n',h0,r);


%% 4.	Hereafter plot the rating curve represented by two straight lines, taking into account that the floodplains start to be inundated at h=4.5m. Therefore the injection point can be selected at h= 4.5 m
% I dont understand the question.
inflex  =   4.5;	% [m]
xdata   = log10([data.Q(data.H>inflex);data.floodQ(data.floodH>inflex)]');
ydata   = log10([data.H(data.H>inflex);data.floodH(data.floodH>inflex)]');

figure(4)
cax     = gca;
set(cax,'NextPlot','add');
xlabel('log10(Q m3/s)');
ylabel('log10(h m)');

% data points scattered
scatter(xdata,ydata);

% Regression line plotted
[r,a,h0_]   = regression(xdata,ydata);
regx        = linspace(0,max(xdata));
regy        = a*regx+h0_;
plot(regx,regy);

h0          = 10^h0_;

fprintf('the h0, the offset is estimated on %f m, which is a lot lower than the not-log\nThe correlation coefficient is %f, which is the same as before\n',h0,r);

%% 5.	What is the equation for the rating curve for each of its branches.

fprintf('The equation for the lineair regression line: y = %f * x + %f\n',a,h0)
